C                                                           PREWIT
C
C  Prewhitens an input sequence in-place.  Uses a low-order prediction error
C    filter. Order selected for stability if needed.
C
C  Author:  George Randall after Dave Harris
C
C  Created: February 12, 1986
C
C  Last Modified:  February 12, 1986
C
C  Input arguments:
C  ----- ----------
C
C    DATA                 REAL*4 array containing input sequence - contains
C                         prewhitened sequence upon exit from this routine.
C
C    #SAMPLES             Number of data points in sequence.
C
C    PREDICTOR_ORDER      Order of prediction filter used to prewhiten
C                         sequence. Truncated for stability
C
C  Output Arguments:
C  ------ ----------
C
C    DATA                 As above.
C
C    A                    Array of prewhitening filter coefficients.
C
C    ERROR_MESSAGE        CHARACTER*130 variable containing error message if
C                         error is detected.  Equal to ' ' if no errors.
C
C  Linkage:     DIRCOR, LEVIN, PEF
C
C
      subroutine prewit(data, nsamps, order, a, errmsg)
C
	parameter ( NCMAX = 12 )
        real*4 data(1), a(1), atemp(NCMAX), sa(NCMAX+1)
        real*4 r(NCMAX), reflct(NCMAX+1)
        character*130 errmsg
        integer order, nsamps, torder
C
C  Initializations
C
        errmsg = ' '
C
C  Range check for predictor order
C
        if (  order .lt. 1 .or. order .gt. NCMAX ) then
          errmsg = '*** PREWIT  -  Predictor order out of bounds (1-12)
     1***'
        end if
C
C  Design the prewhitening filter
C
        call dircor(data, data, nsamps, 0, order+1, r)
        call levin(r, a, reflct, order+1)
C
C  Check for stability, and truncate filter if necessary
C  so the recursive de-whitener is safely stable
C  Use an adaptation of LPTRN code from Markel and Gray
C  found in the IEEE ASSP book of signal processing codes
C
	torder = order
c	do 100 i = NCMAX+1, 1, -1
c	  if ( abs( reflct(i) ) .gt. 0.95 ) torder = i - 1
c 100    continue
C
C  If any reflection coefficients are too close to +1 or -1
C  then the system is getting dangerous, so truncate, and
C  then regenerate the filter coeffs for the truncated filter
C
	if ( torder .ne. order ) then
	  do 110 i = 1, torder
	    sa(i) = reflct(i)
 110        continue
	  do 130 j = 2, torder
	    j2 = j/2
	    q = reflct(j)
	    do 120 k = 1, j2
	      kb = j - k
	      at = sa(k) + q * sa(kb)
	      sa(kb) = sa(kb) + q * sa(k)
	      sa(k) = at
 120          continue
 130        continue
	  do 140 i = 1, torder
	    a(i+1) = sa(i)
 140        continue
	  do 150 i = torder+2, order+1
	    a(i) = 0.
 150	    continue
	  a(1) = 1.
	  order = torder
	endif
C
C  Prewhiten the data
C
C    Negate coefficients of prediction error filter generated by LEVIN.
C    For historical reasons, different storage modes are used in LEVIN
C    and PEF.
C
        do    1 i = 1, order
          atemp(i) = -a(i+1)
C             I
    1   continue
C
        call pef(data, nsamps, atemp, order, 0, 0., data, errmsg)
        if (.not.( errmsg .eq. ' ' )) then
          call append(errmsg, '$ from (PREWIT)$')
        end if
C
      return
      end
